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Abstract 

We present a study of the equilibrium properties of sp-bonded solids within 
the pseudopotential approach, employing recently proposed generalized gradi- 
ent approximation (GGA) exchange correlation functionals. We analyze the 
effects of the gradient corrections on the behavior of the pseudopotentials and 
discuss possible approaches for constructing pseudopotentials self-consistently 
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in the context of gradient corrected functionals. The calculated equilibrium 
properties of solids using the GGA functionals are compared to the ones ob- 
tained through the local density approximation (LDA) and to experimental 
data. A significant improvement over the LDA results is achieved with the 
use of the GGA functionals for cohesive energies. For the lattice constant, the 
same accuracy as in LDA can be obtained when the nonlinear coupling between 
core and valence electrons introduced by the exchange correlation functionals 
is properly taken into account. However, GGA functionals give bulk moduli 
that are too small compared to experiment. 

There have been extensive efforts recently to improve the accuracy of density 
functional theory (DFT) [1] by going beyond the local density approximation (LDA) 
[2]. Including gradient corrections to LDA represents a promising scheme that is 
conceptually simple [3]- [10]. This approach is referred to as Generalized Gradient 
Approximation, or GGA. Calculations using different gradient-corrected functionals 
have been performed to test the applicability of this approach on a variety of systems 
[11]-[17]. 

In an earlier work [17], we investigated one GGA functional recently proposed by 
Perdew and Wang (PW91) [6], [9] in calculations of both atoms and solids. We showed 
that by simply combining the PW91 functional with the pseudopotential approach 
leads to lattice constants for solids, such as simple metals and semiconductors, that 
are larger than experiment, and the percentage errors are significantly larger than 
those obtained from LDA. Results from all-electron calculations [13], [16] for systems 
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other than the ones we considered, such as transition metals, suggest that there is no 
such significant increase of the error in lattice constants when the PW91 functional 
is used. In order to understand this difference between the results obtained from 
the pseudopotential and the all-electron calculations, we first examined the charge 
distribution of individual atoms. We considered both PW91 and the more recently 
proposed functional by Lacks and Gordon (LG) [10]. The reasons for choosing these 
two functionals are: (1) They are among the most recent additions to the list of pro- 
posed gradient corrections and are intended to give better results compared to earlier 
attempts. (2) The two functionals have similar expressions which simplifies compu- 
tational implementation. (3) The Lacks-Gordon exchange functional is produced by 
fitting to exact results, while the PW91 functional is derived from first principles; 
comparison between these two functionals may provide insight for a more accurate 
approach. 

In studying the charge distributions of individual atoms, we have found it instruc- 
tive to evaluate the following quantity 



C n i(R) is the difference between the charge enclosed within a sphere of radius R 
around the nucleus calculated from the LDA and the GGA functional for each single 
electron orbital nl. In Fig. 1 we display C n i(R) for the nl = 2p, 3s and 3p orbitals 
of Si, and the nl = 2p, 3s, 3p orbitals of Na. The more positive C n i{R) is the more 
charge has been pushed outside the region contained by the sphere of radius R in the 
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GGA calculation compared to the charge obtained by LDA. 

In order to see the difference in the physically important range, we used the 
experimentally measured bond length as the unit along the x-axis in Fig. 1. As 
is obvious from this figure, there is almost no difference in the charge distribution 
between LDA and GGA results for the 2p core orbital of Si. For the 3s and 3p 
valence orbitals of Si, in the neighborhood of the bonding region substantial charge 
has been pushed away from the nucleus in both the PW91 and the LG calculations, 
relative to the LDA results. A similar situation occurs in the case of Na as is shown 
in Fig. 1. 

These comparisons indicate a weaker interaction between the valence electrons 
of the atom and the ion when using GGA functionals as opposed to LDA, due to 
the spreading out of the valence charge in GGA calculations. In an approach which 
simply replaces the effects due to the ion and the core electrons with a pseudo- 
core which is constructed to reproduce the results of the all-electron calculation, 
this will unavoidably lead to a weaker interaction between the valence electrons and 
the pseudo-core. Since in the pseudopotential framework the properties of solids 
are determined by the interaction between valence electrons and the pseudo-core, 
the tendency for valence charge to be pushed away from the nucleus when GGA 
functionals are used leads to a softer solid, characterized by larger equilibrium lattice 
constant and smaller bulk modulus than LDA. 

This observation points to the necessity of properly taking into account the non- 



4 



linear coupling between the valence and core electrons in the exchange correlation 
functional within the pseudopotential approach, when GGA functional are used. To 
examine the effects due to the inclusion of gradient corrections, we first consider the 
behavior of the quantity s which is used in the definition of GGA functionals in ad- 
dition to the charge density n. s is defined as a function of the charge density n and 
its gradient Vn: 



Although the charge density n can be separated into the core charge density n c and 
the valence charge density n v so that n = n c + n v , such a separation can not be 
written for s due to the nonlinearity of the expression of Eq. (2). We use Si as an 
example to further illustrate this point. We display in Fig. 2 both the core and valence 
charge density for the Si atom, calculated with the LDA and the GGA functionals 
respectively. There is evidently no significant difference between the results of these 
two calculations as far as the overlap between charge density is concerned. Since the 
results from the PW91 and LG functionals tend to be very similar (see for example 
the comparison in Fig. 1), from now on we will only present results obtained with 
PW91. Fig. 3(a) displays the quantity 



which is a measure of the nonlinearity of the exchange- correlation potential. There 
is a substantial increase in the nonlinearity of the exchange-correlation potential in 
the region where the overlap between the core and valence charge is not negligible 




(2) 



A\4 C = V xc [n c ] + V xc [n v ] - V xc [(n c + n v )\ 



(3) 
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(compare with Fig. 2). This increase in nonlinearity is due to the fact that the 
variables used in the expression for the GGA functionals are not separable in terms 
of the valence and core parts. For example, we show in Fig. 3(b) the values of s c , 
s v and s which correspond to the value of the quantity defined in Eq. (2) calculated 
with the core, valence, and total charge density respectively. Therefore, the inclusion 
of gradient corrections seems to increase the effects of coupling between the valence 
and core electrons. 

As was discussed in the previous paragraphs, the simple unscreening of the pseu- 
dopotential by 

Kon = V screened ~ VffM - V xc [n v ] (4) 

where V# and V xc are the Hartree and the exchange-correlation potentials respectively, 
does not give satisfactory results for the properties of solids. This approximation 
implicitly assumes the linearization of the exchange-correlation functional. As was 
pointed out by Louie el al. [18], a more consistent approach is to include the core 
charge density in the unscreening. That is, instead of taking out only the \4 c [n„] part 
in the unscreening procedure as in Eq. (4), V^ s n should be defined instead as: 

Vion = ^screened ~ V»M - V xc [(n c + U v )} (5) 

where n c is a rigid core charge density, constructed from a reference atomic system. 
Since the core electrons are not included in the pseudopotential calculation, when- 
ever the exchange-correlation energy and potential are needed, the full charge density 
n — n c + n v must be used. This procedure is exact within the rigid core approxima- 
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tion, but it would require a very large number of plane waves to describe the core 
charge density accurately, and one loses the advantages of using the pseudopotential 
formalism by adopting this approach. So even though it is theoretically correct, it is 
not practical from a computational point of view. Therefore, it is necessary to make 
some approximation in order to obtain a practical computational scheme. 

In the present paper, we follow the partial core prescription proposed in Ref. [18]. 
The full core charge density is replaced with an artificial core charge density h c defined 
as: 

A sin Br 

n c = ,r<r c 

r 

h c = n c ,r>r c (6) 

The parameters A and B are determined by the requirement that the value of n c 
and its derivative with respect to the radius r be exactly the same as those of the 
real core charge density n c at the cutoff radius r c . We have found that in order to 
capture the nonlinear coupling between the core and the valence electrons for the case 
of the PW91 functional, it is necessary to use r c smaller than what was suggested 
in Ref. [18] for LDA calculations. In our calculations, r c is chosen as the radius 
where the core charge density n c is 6-7 times larger than the valence charge density 
n v . In Ref. [18] (which dealt with LDA calculations), r c was chosen as the radius 
where the core charge density n c was 2-3 times larger than the valence charge density. 
It is worthwhile mentioning that the pathological oscillatory behavior of the PW91 
exchange-correlation potential near the nuclei, which causes problems in creating 
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smooth pseudopotentials during the unscreening procedure [13], [17], is automatically 
eliminated by using the partial core correction. 

As an illustration of how this approach works, we consider four sp-bonded solids 
in their ground state phase: Si (diamond), Ge (diamond), GaAs (zincblende), and 
Al (fee). We construct self-consistent pseudopotentials as described in Ref. [17] with 
the partial core prescription for the unscreening procedure discussed above. For the 
LDA calculations we use the exchange-correlation potential of Ceperley and Alder 
as parametrized by Perdew and Zunger [19] and norm-conserving pseudopotentials 
from Bachelet, Hamann, and Schluter (BHS) [20]. We use a plane-wave basis for 
the expansion of the wavefunction of valence electrons: the highest kinetic energy of 
the plane waves in the basis is 16 Ry. For reciprocal space integrations, 29 special 
k-points in the irreducible Brillouin zone are used for the diamond structure and 
the zinc-blende structure, and 213 special k-points for the fee structure. The gradient 
and the laplacian of the density, which are needed for the GGA functionals considered 
here, were obtained through FFT's, with minimal increase in CPU time (less than 
3%). The calculated energy vs. volume results are fitted to the Birch- Murnaghan 
equation of state [21]. The equilibrium properties are then derived from the equation 
of state curves. The cohesive energy is taken to be the total energy difference between 
the solid in equilibrium and the isolated atom. Spin polarization effects on the free 
atom energy are taken into account with the empirical formula AE P = —0.18 x rip, 
[22], [23], where n p = n-\ — with the number of electrons having spin up and 
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n-i the number of electrons with spin down. Spin polarization is expected to have 
negligible effects on the total energy E of nonmagnetic solids. 

The calculated ground state properties using LDA and PW91 are summarized 
and compared to experimental data in Table I. The results from PW91 represent a 
substantial improvement of the over-binding problem of LDA: the cohesive energies 
are in better agreement with experiment for all the solids we have considered. This 
improvement in cohesive energies can be attributed to a large extene to the fact 
PW91 gives a more accurate atomic energy than LDA [17]. For the equilibrium lattice 
constant, the value obtained from PW91 is consistently larger than LDA results. In 
the case of Al, this makes the value obtained from PW91 closer to experiment than 
the LDA result. For Si, Ge, and GaAs, the results obtained from LDA and PW91 
are of the same accuracy compared to experiment. PW91 tends to overcorrect the 
LDA results and gives an overestimate for the equilibrium lattice constant of these 
systems. For the bulk modulus, the values obtained from PW91 are smaller than 
the LDA results. While this leads to a better result for Al, the bulk moduli we 
obtained for the three semiconductors are significantly underestimated (by -12 % to 
-25 % compared to experiment). Similar observations for earlier gradient-corrected 
functionals have been reported [11], [12]. Finally, we compare our results to recent all- 
electron, linearized augamented planewave (LAPW) [24] total energy and electronic 
structure calculations, with the same GGA functional as in our work. It is obvious 
from the comparison of Table I that the present pseudopotential calculation results 
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with the partial core correction represent a significant improvement over the results 
without the partial core correction, and agree very well with the all-electron LAPW 
calculation results. The remaining discrepancy is probably due to the relaxation of 
the core electrons which is not allowed in the pseudopotential calculation. 

For the electronic structure, we compare in Table II the band gap predicted by 
LDA and PW91 at both the experimentally measured lattice constant and the theo- 
retical equilibrium lattice constant. An overall improvement which brings the values 
closer to the experimental data at the experimental lattice constant is found for all 
three semiconductors we have considered, although the magnitude of the improvement 
depends on the material. There is no consistent improvement for the band gaps at the 
theoretical equilibrium lattice constant. We note here that this is simply a comparison 
between LDA and PW91 as different approximations to the exchange correlation func- 
tionals. The well known inability of density functional theory to reproduce accurately 
band gaps in semiconductors and insulators is much more complicated and related 
to the intrinsic discontinuity of the exchange correlation functional [25] , which is not 
represented by either of the two approximations used here. Good agreement for the 
band gap values nevertheless can be obtained by using DFT/LDA wavefunctions and 
solving the self-energy operator equations, within the so called GW approximation 
[26]. 

In conclusion, we showed that it is essential to take into account the core-valence 
coupling in the pseudopotential calculations when using GGA exchange-correlation 
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functionals. To this end, we have found that the partial core prescription of Louie et 
al. [18] is most appropriate when using a plane-wave basis. We considered the struc- 
tural properties of Si, Ge, GaAs, and Al using both LDA and PW91. We found that 
PW91 gives consistently better cohesive energies than LDA. We also demonstrated 
that for the lattice constant the same accuracy as in LDA can be obtained with GGA, 
as long as the nonlinearity of the gradient corrected functional is properly taken into 
account. For the semiconductors we considered, the bulk moduli obtained with the use 
of GGA functionals represent significant underestimates of the experimental results. 
The PW91 functional does give a better description for the equilibrium properties 
of Al. We conclude that further search may be needed for an exchange-correlation 
functional which is consistently better for all solids. 

In view of the above results, one may inquire what are the physical situations in 
which the use of GGA functionals can provide significant improvements over LDA 
results. Recently, calculations have been reported for H 2 dissociation on a Cu(lll) 
surface with the LDA and the GGA [27], [28]. The GGA results for this system rep- 
resent significant improvements over the LDA results. It has also been demonstrated 
that the GGA gives results in better agreement with experiments than the LDA for 
finite systems (atoms and molecules) and metallic surfaces [13]- [15]. It is therefore 
expected that the GGA will give, in general, a better description for the interaction 
between molecules and other molecules or solid surfaces. The reason that the GGA 
should give better results for these interactions can be attributed to the fact that 
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substantial part of the interactions in the these systems are related to the tails of the 
electronic wave functions, where the GGA gives a more accurate description than the 
LDA. 

This work was supported by the Materials Research Laboratory of Harvard Uni- 
versity which is funded by National Science Foundation Grant # DMR 89-20490. The 
calculations were carried out at the Cornell National Supercomputer Facility. 
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Figure Captions 

1. Integrated charge difference C n i(R) [see Eq. (1)] calculated with PW91 (solid 
lines), and LG (dashed lines) for nl =2p, 3s, 3p orbitals of Si, and nl =2p, 3s, 
3p orbitals of Na. The results for the 3p orbital of Na have been divided by a 
factor 2 so that they can be displayed on the same scale as the results for Si 

2. The core (upper panel) and valence (lower panel) charge density of the silicon 
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atom calculated with LDA and GGA functional as a function of distance from 
the nucleus. Notice the different density scales in the two panels. 

3. (a) The quantity AV XC [see Eq. (3)] for the silicon atom calculated with LDA 
and PW91 as a function of distance from the nucleus. 

(b)The values of s [see Eq. [2]] calculated from the core (s c ), the valence (s„), 
and the total charge density (s) for the silicon atom as a function of distance 
from the nucleus. 
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